/*
    Copyright 2006-2008 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <iostream.h>
#include <math.h>
#include "grid.h"
#include "optical.h"
#include "sphgrid.h"
#include "sphparticle.h"
#include "misc.h"
#include "counter.h"
#include "integrator2.h"

class sphintegrator : public integrator2 {
public:
  class yintegrator : public integrator2 {
  public:
    class zintegrator : public integrator2 {
    public:
      sphparticle* p;
      double x,y;
      double integrand(double); // The integrand
    };
    double integrand(double);
    zintegrator zi;
    double zlo,zhi;
  };
  double integrand(double);
  yintegrator yi;
  double ylo,yhi;
  double project(sphparticle*,grid_cell<optical_data>*);
};

double sphintegrator::project(sphparticle* p, 
			      grid_cell<optical_data>* c)
{
  // Avoid integrating where the kernel is identically zero
  // so check what the overlapping cube is
  double pminx=p->pos[0]-2*p->h;
  double pminy=p->pos[1]-2*p->h;
  double pminz=p->pos[2]-2*p->h;

  double pmaxx=p->pos[0]+2*p->h;
  double pmaxy=p->pos[1]+2*p->h;
  double pmaxz=p->pos[2]+2*p->h;

  vec3<gt> cellmin=c->getmin();
  vec3<gt> cellmax=c->getmax();
  int notwithin=0;
    
  if(cellmin[0]>pminx) {
    pminx=cellmin[0];
    notwithin=1;
  }
  if(cellmin[1]>pminy) {
    pminy=cellmin[1];
    notwithin=1;
  }
  if(cellmin[2]>pminz) {
    pminz=cellmin[2];
    notwithin=1;
  }
  
  if(cellmax[0]<pmaxx) {
    pmaxx=cellmax[0];
    notwithin=1;
  }
  if(cellmax[1]<pmaxy) {
    pmaxy=cellmax[1];
    notwithin=1;
  }
  if(cellmax[2]<pmaxz) {
    pmaxz=cellmax[2];
    notwithin=1;
  }
  
  if( notwithin && (pmaxx>pminx) && (pmaxy>pminy) && (pmaxz>pminz) ) {

    // We have a nonzero region where kernel and cell overlap - integrate
    ylo=pminy;
    yhi=pmaxy;
    yi.zlo=pminz;
    yi.zhi=pmaxz;
    yi.zi.p=p;

    // Because of how the integrator works we don't want the identically
    // zero region in the beginning since then the step size is increased and
    // overshoots when the kernel goes nonzero leading to a huge penalty
    // so, we divide the integral in two 
    if(pminx>= p->pos[0])
      return integrate(pminx,pmaxx);
    else if(pmaxx<= p->pos[0])
      return -integrate(pmaxx,pminx);
    else
      return -integrate((double)p->pos[0],pminx)+
	integrate((double)p->pos[0],pmaxx);
    
  }
  else if( notwithin )
    return 0;
  else
    return 1;
    
}

double sphintegrator::integrand(double x)
{
  yi.zi.x=x;

    if(ylo>=yi.zi.p->pos[1])
      return yi.integrate(ylo,yhi);
    else if(yhi<=yi.zi.p->pos[1])
      return -yi.integrate(yhi,ylo);
    else
      return -yi.integrate((double)yi.zi.p->pos[1],ylo)+
	yi.integrate((double)yi.zi.p->pos[1],yhi);
  
}

double sphintegrator::yintegrator::integrand(double y)
{
  zi.y=y;


    if(zlo>=zi.p->pos[2])
      return zi.integrate(zlo,zhi);
    else if(zhi<=zi.p->pos[2])
      return -zi.integrate(zhi,zlo);
    else
      return -zi.integrate((double)zi.p->pos[2],zlo)+
	zi.integrate((double)zi.p->pos[2],zhi);
}

double sphintegrator::yintegrator::zintegrator::integrand(double z)
{
  double px=p->pos[0];
  double py=p->pos[1];
  double pz=p->pos[2];

  // Here's a function of the coordinates x,y,z and the sphparticle p
  double r=sqrt( (x-px)*(x-px) + (y-py)*(y-py) + (z-pz)*(z-pz) );

  return p->kernel(r);
}

void makeinterpolationdata();
void testinterpolationdata();
  sphintegrator sphint;

main()
{
  sphint.setaccuracy(1e-4,1e-4);
  sphint.yi.setaccuracy(1e-4,1e-4);
  sphint.yi.zi.setaccuracy(1e-4,1e-4);

  //makeinterpolationdata();
  testinterpolationdata();
}

void testinterpolationdata()
{
  ifstream inf("/home/sunrise/dust/SNAP010/SNAP010-550.sph");
  double q;
  inf >> q >> q >> q;
  sphparticle p;
  p.h=0.25;
  p.pos=vec3<double>(0.5,0,0);
  sphgrid g(-1.5,4.5,2,-1.5,1.5,1,-1.5,1.5,1);
  grid_cell<optical_data>* c =g.locate(vec3<double>(0,0,0));
  g.refine(c,vec3<int>(2,2,2));
  //c =g.locate(vec3<double>(0,0,0));
  //->refine(c,vec3<int>(2,2,2));

  while(inf.good()) {
    inf >> p;
    g.init();
    while(c=g.enumerate()) {
      double a1=p.project2(c);
      double a2=sphint.project(&p,c);
      double r=(c->getcenter()-p.pos).mag();
      double a3=p.kernel(r)*c->volume();
      if(a2!=0)
	cout << a1 << '\t' << a2 << '\t' << a3 << '\t' << p.h << '\t' 
	     << r << '\t' << c->getsize()[0] << '\n';
    }
    }
}

void makeinterpolationdata()
{
  int nx=13;
  int ny=13;
  int nz=13;
  int nh=20;

  cout << "4\n" << nx << ' ' << ny << ' ' << nz << ' ' << nh << "\n";

  double xx[nx];
  double yy[ny];
  double zz[nz];
  double hh[nh];

  int i=1;
  xx[0]=-1e30;
  xx[nx-1]=1e30;
  for(double x=-2;x<=2.01;x+=4./(nx-3)) {
    xx[i++]=x;
  }

  i=1;
  yy[0]=-1e30;
  yy[ny-1]=1e30;
  for(double y=-2;y<=2.01;y+=4./(ny-3)) {
    yy[i++]=y;
  }

  i=1;
  zz[0]=-1e30;
  zz[nz-1]=1e30;
  for(double z=-2;z<=2.01;z+=4./(nz-3)) {
    zz[i++]=z;
  }

  i=1;
  hh[0]=0.0001;
  for(int l=1;l<nh;l++) {    
    hh[l]=0.01*pow(10.,4.*(l-1)/(nh-2));
  }

  for(int i=0;i<nx;i++)
    cout << xx[i] << ' ';
  cout << '\n';
  for(int j=0;j<ny;j++)
    cout << yy[j] << ' ';
  cout << '\n';
  for(int k=0;k<nz;k++)
    cout << zz[k] << ' ';
  cout << '\n';
  for(int l=0;l<nh;l++)
    cout << hh[l] << ' ';
  cout << "\n\n";
    
  sphgrid g(-0.5,0.5,1,-0.5,0.5,1,-0.5,0.5,1);


  sphparticle p;
  p.m_g=1;
  p.h=0.0001;
  p.pos=vec3<gt>(0,0.5,0);
  grid_cell<optical_data>* c =g.locate(vec3<double>(0.));
  //cout << sphint.project(&p,c) << '\n';

  double y=0;
  counter cc(10);
  for(int l=7;l<nh;l++)
    for(int k=1;k<nz;k++)
      for(int j=1;j<ny;j++) {
	for(int i=6;i<nx;i++) {
	  //cerr << l << ' ' << k << ' ' << j << "    \r";
	  //cc.incr();
	  p.h=hh[l];
	  p.pos[0]=xx[i]*p.h+0.5;
	  p.pos[1]=yy[j]*p.h+0.5;
	  p.pos[2]=zz[k]*p.h+0.5;
	  if(p.pos[0]<0)
	    p.pos[0]=0;
	  if(p.pos[1]<0)
	    p.pos[1]=0;
	  if(p.pos[2]<0)
	    p.pos[2]=0;
	  double a;
	  a = sphint.project(&p,c);
	  cout << a << ' ';
	  cerr << xx[i] << ' ' << yy[j] << ' ' << zz[k] << ' ' <<p.pos << '\t' << p.h << '\n';
	}
	cout << '\n';
      }
  /*
 
  grid<double> g(-5,5,10,-5,5,10,-5,5,10);
  grid_cell<double>* c;
  c=g.locate(matrix<double>(0.5,0.5,0.5));
  c->refine(matrix<int>(10.,10.,10.));
  c=g.locate(matrix<double>(-.5,0.5,0.5));
  c->refine(matrix<int>(10.,10.,10.));
  
  c=g.locate(matrix<double>(0.05,0.55,0.55));
  c->refine(matrix<int>(10.,10.,10.));
  c=g.locate(matrix<double>(-0.05,0.55,0.55));
  c->refine(matrix<int>(10.,10.,10.));
  
  matrix<double> p(0.005,0.555,0.555);

  grid_cell<double>* d;
  grid_cell<double>* dd;
  int t1,t2,t3,t4;
  for(int i=0;i<10;i++) {
    c=g.locate(p);
    t1=realcc();
    d=c->neighbour(p,0,1);
    t2=realcc();
    dd=c->neighbour(p,0,-1);
    t3=realcc();
    t4=realcc();
    }
  cout << t2-t1 << '\t' << t3-t2 << '\n';
  cout << c->getmin() << c->getmax() << '\n';
  cout << d->getmin() << d->getmax() << '\n';
  cout << dd->getmin() << dd->getmax() << '\n';
  */  
  /*  
  const int n=300;
  cout << n+1 << '\n';
  for(double i=0;i<=1;i+=1./n)
    cout << i << '\t' << cos(2*pi*i) << '\n';
  */

  //linear_interpolator i("sin.inp");
  //cout << i.interpolate(0.654);
}
